*Set root data directory
local rootdir
cd "`rootdir'"

set maxvar 12000
set more off
set matsize 11000
set emptycells drop


foreach round in 38{
foreach i in grains veg{


use round`round'_edit, clear
merge 1:1 sector subround region subsample segment substratum fsu hhno using round`round'_hh, keepusing(mpce hhsize mult d*)
drop _merge
merge m:m region using mult`round', keepusing(regioncode)
drop _merge

egen hhid=group(sector subround region subsample segment substratum fsu hhno)


drop hq* hx* *intox* *cloth* *light*

gen statecode=floor(region/10)

sort sector region fsu

merge m:m statecode using nss50_nm61
drop if fsu==.
drop if _merge==2
drop _merge

ren statename statenamedist

drop if fsu==.
drop if statenamedist==""


keep *`i'* subround mpce round regioncode sector statenamedist fsu hhid hhsize

keep if sector==1

gen fn=0
gen rn=0
gen sn=0


levelsof n`i'tot, local(nmax)

forvalues j=1(1)`nmax'{
gen uv`i'_`j'=x`i'_`j'/q`i'_`j'
replace uv`i'_`j'=. if uv`i'_`j'<=0
egen p99`i'_`j'=pctile(uv`i'_`j'), p(99)
egen p1`i'_`j'=pctile(uv`i'_`j'), p(1)
replace uv`i'_`j'=p99`i'_`j' if uv`i'_`j'>p99`i'_`j' & uv`i'_`j'~=.
replace uv`i'_`j'=p1`i'_`j' if uv`i'_`j'<p1`i'_`j' & uv`i'_`j'~=.
drop p99* p1*

bysort statenamedist fsu: egen fmeduv`i'_`j'=median(uv`i'_`j')
bysort regioncode: egen rmeduv`i'_`j'=median(uv`i'_`j')
bysort statenamedist: egen smeduv`i'_`j'=median(uv`i'_`j')


bysort statenamedist subround: egen sr_smeduv`i'_`j'=median(uv`i'_`j')
bysort regioncode: egen sr_rmeduv`i'_`j'=median(uv`i'_`j')

gen hpresence`i'_`j'=0
replace hpresence`i'_`j'=1 if x`i'_`j'>0 & fmeduv`i'_`j'~=.
gen fpresence`i'_`j'=0
replace fpresence`i'_`j'=1 if fmeduv`i'_`j'~=.
gen rpresence`i'_`j'=0
replace rpresence`i'_`j'=1 if rmeduv`i'_`j'~=.
gen spresence`i'_`j'=0
replace spresence`i'_`j'=1 if smeduv`i'_`j'~=.

replace fn=fn+fpresence`i'_`j'
replace rn=rn+rpresence`i'_`j'
replace sn=sn+spresence`i'_`j'

bysort statenamedist fsu: egen fnum`i'_`j'=total(fpresence`i'_`j')
bysort regioncode: egen rnum`i'_`j'=total(rpresence`i'_`j')
bysort statenamedist: egen snum`i'_`j'=total(spresence`i'_`j')

bysort statenamedist fsu: egen fx`i'_`j'=total(x`i'_`j')
bysort regioncode: egen rx`i'_`j'=total(x`i'_`j')
bysort statenamedist: egen sx`i'_`j'=total(x`i'_`j')
drop uv`i'_`j' 
}

bysort statenamedist fsu: gen fnumhh=_N
bysort regioncode: gen rnumhh=_N
bysort statenamedist: gen snumhh=_N

bysort statenamedist fsu: egen ftotx=total(x`i')
bysort regioncode: egen rtotx=total(x`i')
bysort statenamedist: egen stotx=total(x`i')

bysort statenamedist fsu: egen fmaxx=max(x`i')
bysort regioncode: egen rmaxx=max(x`i')
bysort statenamedist: egen smaxx=max(x`i')

duplicates drop hhid, force

expand `nmax'
bysort hhid: gen gcounter=_n
foreach k in sr_smeduv sr_rmeduv fmeduv rmeduv smeduv hpresence fpresence rpresence spresence fnum rnum snum x q fx rx sx{
gen `k'=.
forvalues j=1(1)`nmax'{
replace `k'=`k'`i'_`j' if gcounter==`j'
drop `k'`i'_`j'
}
}



save temp`i'`round', replace
}
}




foreach round in 66{
foreach i in grains veg{


use round`round'_edit, clear
merge 1:1 sector subround region subsample segment substratum fsu hhno using round`round'_hh, keepusing(mpce hhsize mult d*)
drop _merge
merge m:m region using mult`round', keepusing(regioncode)
drop _merge

egen hhid=group(sector subround region subsample segment substratum fsu hhno)

drop hq* hx* *intox* *cloth* *light*
gen statecode=floor(region/10)
gen districtcode=district

sort sector region fsu
bysort fsu sector: egen minddist=min(districtcode)
bysort fsu sector: egen maxddist=max(districtcode)
replace districtcode=minddist if districtcode==.
replace districtcode=maxddist if districtcode==0 & maxddist~=.

merge m:m statecode districtcode using nss`round'_nm61"drop if fsu==.
drop if _merge==2
drop _merge

ren statename statename`round'
ren districtname districtname`round'

merge m:m statename`round' districtname`round' using district_conversion
drop if fsu==.
drop if _merge==2
drop _merge

drop if statenamedist==""
drop if districtnamedist<"Adilabad"


keep *`i'* subround region regioncode mpce round sector statenamedist districtnamedist fsu hhid hhsize

gen fn=0
gen rn=0
gen dn=0
gen sn=0


keep if sector==1

levelsof n`i'tot, local(nmax)

forvalues j=1(1)`nmax'{
gen uv`i'_`j'=x`i'_`j'/q`i'_`j'
replace uv`i'_`j'=. if uv`i'_`j'<=0
egen p99`i'_`j'=pctile(uv`i'_`j'), p(99)
egen p1`i'_`j'=pctile(uv`i'_`j'), p(1)
replace uv`i'_`j'=p99`i'_`j' if uv`i'_`j'>p99`i'_`j' & uv`i'_`j'~=.
replace uv`i'_`j'=p1`i'_`j' if uv`i'_`j'<p1`i'_`j' & uv`i'_`j'~=.
drop p99* p1*

*unit value is also a measure of availability here
bysort statenamedist districtnamedist fsu: egen fmeduv`i'_`j'=median(uv`i'_`j')
bysort statenamedist districtnamedist: egen dmeduv`i'_`j'=median(uv`i'_`j')
bysort regioncode: egen rmeduv`i'_`j'=median(uv`i'_`j')
bysort statenamedist: egen smeduv`i'_`j'=median(uv`i'_`j')

bysort statenamedist districtnamedist subround: egen sr_dmeduv`i'_`j'=median(uv`i'_`j')
bysort regioncode: egen sr_rmeduv`i'_`j'=median(uv`i'_`j')
bysort statenamedist subround: egen sr_smeduv`i'_`j'=median(uv`i'_`j')

gen hpresence`i'_`j'=0
replace hpresence`i'_`j'=1 if x`i'_`j'>0 & fmeduv`i'_`j'~=.
gen fpresence`i'_`j'=0
replace fpresence`i'_`j'=1 if fmeduv`i'_`j'~=.
gen dpresence`i'_`j'=0
replace dpresence`i'_`j'=1 if dmeduv`i'_`j'~=.
gen rpresence`i'_`j'=0
replace rpresence`i'_`j'=1 if rmeduv`i'_`j'~=.
gen spresence`i'_`j'=0
replace spresence`i'_`j'=1 if smeduv`i'_`j'~=.

replace fn=fn+fpresence`i'_`j'
replace dn=dn+dpresence`i'_`j'
replace rn=rn+rpresence`i'_`j'
replace sn=sn+spresence`i'_`j'

bysort statenamedist districtnamedist fsu: egen fnum`i'_`j'=total(fpresence`i'_`j')
bysort statenamedist districtnamedist: egen dnum`i'_`j'=total(dpresence`i'_`j')
bysort regioncode: egen rnum`i'_`j'=total(rpresence`i'_`j')
bysort statenamedist: egen snum`i'_`j'=total(spresence`i'_`j')

bysort statenamedist districtnamedist fsu: egen fx`i'_`j'=total(x`i'_`j')
bysort statenamedist districtnamedist: egen dx`i'_`j'=total(x`i'_`j')
bysort regioncode: egen rx`i'_`j'=total(x`i'_`j')
bysort statenamedist: egen sx`i'_`j'=total(x`i'_`j')
drop uv`i'_`j'
* q`i'_`j'
}

bysort statenamedist districtnamedist fsu: gen fnumhh=_N
bysort statenamedist districtnamedist: gen dnumhh=_N
bysort regioncode: gen rnumhh=_N
bysort statenamedist: gen snumhh=_N

bysort statenamedist districtnamedist fsu: egen ftotx=total(x`i')
bysort statenamedist districtnamedist: egen dtotx=total(x`i')
bysort regioncode: egen rtotx=total(x`i')
bysort statenamedist: egen stotx=total(x`i')

bysort statenamedist districtnamedist fsu: egen fmaxx=max(x`i')
bysort statenamedist districtnamedist: egen dmaxx=max(x`i')
bysort regioncode: egen rmaxx=max(x`i')
bysort statenamedist: egen smaxx=max(x`i')

duplicates drop hhid, force

expand `nmax'
bysort hhid: gen gcounter=_n
foreach k in sr_dmeduv sr_rmeduv sr_smeduv fmeduv dmeduv rmeduv smeduv hpresence fpresence dpresence rpresence spresence fnum dnum rnum snum x q fx dx rx sx{
gen `k'=.
forvalues j=1(1)`nmax'{
replace `k'=`k'`i'_`j' if gcounter==`j'
drop `k'`i'_`j'
}
}



save temp`i'`round', replace
}
}





*****Elasticity**********




***Region (all rounds)

foreach i in grains veg{

local round

use temp`i'38, clear
append using temp`i'66


gen year=1983 if round==38
replace year=2009 if round==66

**merge in rainfall
merge m:1 regioncode year using "D:\Dropbox\indiadistrict\region_rain"
gen rain=act_rain_gca-norm_rain_gca

**region drop for rainfall analysis
keep if rain~=.
keep if _merge==3
drop if regioncode==.

*rural only
keep if sector==1
drop if reglat==. | reglon==.


bysort regioncode round: gen counta=_n
replace counta=0 if counta~=1
by regioncode: egen totcounta=total(counta)
drop counta


*threshold for share of housheolds purchasing
drop rnumhh
bysort round regioncode gcounter: egen frac_hh=total(hpresence)
bysort round regioncode gcounter: gen rnumhh=_N
replace frac_hh=frac_hh/rnumhh


drop rx rtotx
bysort round regioncode gcounter: egen rx=total(x)
bysort round regioncode: egen rtotx=total(x)
gen rshare=rx/rtotx

**drop households with missing data/implausible values
egen hh=group(round regioncode fsu hhid)

drop if x`i'==0
drop if hhsize>14
drop if mpce>10000 | mpce<10
drop if x==0
gen loghhn=log(n`i')

**only keep households with at least 2 varieties
drop if loghhn==0 | loghhn==.


compress



gen share=x/x`i'

replace mpce=mpce/7.256 if round==66
egen medpce=median(mpce)

keep n`i'tot round* regioncode rn gcounter hpresence share rmeduv x`i' hhsize rain hh frac_hh reglat reglon mpce totcounta



egen unit_all=group(gcounter regioncode)
bysort round regioncode gcounter: gen rcounter=_n


levelsof n`i'tot, local(nmax)



gen state=floor(regioncode/10)
gen lns=log(share)


**make base the most popular variety in each region,

bysort regioncode gcounter: egen totx=max(frac_hh)
by regioncode: egen maxx=max(totx)
gen gmax=gcounter if totx==maxx
by regioncode: egen gbase=min(gmax)
summ gbase


gen blns=lns if gcounter==gbase
bysort hh: egen base_lns=min(blns)
gen brmeduv=rmeduv if gcounter==gbase
bysort hh: egen base_rmeduv=min(brmeduv)
gen bfrac_hh=frac_hh if gcounter==gbase
bysort hh: egen base_frac_hh=min(bfrac_hh)


gen dlns=lns-base_lns if lns~=0 & base_lns~=0 & gcounter~=gbase
gen dlnp=log(rmeduv)-log(base_rmeduv) if lns~=0 & base_lns~=0 & gcounter~=gbase
gen dfrac_hh=frac_hh-base_frac_hh if lns~=0 & base_lns~=0 & gcounter~=gbase




cap: drop base_* 
cap: drop blns 
cap: drop brmeduv 
cap: drop gmax
cap: drop maxx
cap: drop totx
cap: drop blnrs 

gen lnxpc=log(x`i'/(hhsize))
gen lnhhsize=log(hhsize)

foreach j in rmeduv{
replace `j'=log(`j')
}

gen prmeduv=rmeduv



foreach j in 66{
gen round`j'=0
replace round`j'=1 if round==`j'
}


egen region=group(regioncode)
egen regionmax=max(region)



sort round regioncode 
local regmax=regionmax[1]
forvalues j=1(1)`regmax'{
gen raintemp=rain if region==`j'
gen lat=reglat if region==`j'
gen lon=reglon if region==`j'
egen tlat=min(lat)
egen tlon=min(lon)
geodist reglat reglon tlat tlon, gen(dist`j') sphere
by round: egen rain_`j'=min(raintemp)
by round regioncode: egen distrain_`j'=min(rain_`j'/(1+log(1+dist`j')))
drop raintemp dist`j' lat lon tlat tlon rain_`j'
}



gen logrn=log(rn)


forvalues j=1(1)`nmax'{
gen pres=1 if gcounter==`j' & hpresence==1
bysort hh: egen pres`j'=min(pres)
drop pres
replace pres`j'=0 if pres`j'==.

gen price=prmeduv if gcounter==`j'
bysort round regioncode: egen price`j'=min(price)
bysort round state: egen price`j'_s=max(price)
bysort round: egen price`j'_ai=max(price)
replace price`j'=price`j'_s if price`j'==.
replace price`j'=price`j'_ai if price`j'==.
drop price price`j'_s price`j'_ai
}






egen unit_ind=group(regioncode)

xtset unit_all


cap: drop totcount
gen count=0
replace count=1 if lns~=.
bysort unit_all count round: gen r2counter=_n
replace r2counter=0 if r2counter~=1
bysort unit_all count: egen totcount=total(r2counter)


gen select=0
replace select=1 if dlns~=.


gen dfrac_hh2=dfrac_hh^2
gen dfrac_hh3=dfrac_hh^3
gen frac_hh2=frac_hh^2
gen frac_hh3=frac_hh^3


cap: drop zzz*
cap: drop gdum*
cap: drop g2dum*


drop if dlns==.


**construct a.good * rain dummies  AND   b. state * good *  rain dummies

egen zzz=group(gcounter)  

egen zzzmax=max(zzz)
local zzmax=zzzmax[1]
forvalues j=2(1)`zzmax'{
gen gdum`j'=0
replace gdum`j'=1 if zzz==`j'
gen g2dum`j'_rain=gdum`j'*rain
}
drop zzz zzzmax gdum*

egen zzz=group(gcounter state)  
egen zzzmax=max(zzz)
local zzmax=zzzmax[1]
forvalues j=2(1)`zzmax'{
gen gdum`j'=0
replace gdum`j'=1 if zzz==`j'
gen g3dum`j'_rain=gdum`j'*rain
}


bysort round regioncode gcounter: egen dlns_mean=mean(dlns)
by round regioncode gcounter: gen mean_counter=_n
by round regioncode gcounter: gen weight_counter=_N
egen unit=group(round regioncode gcounter)

keep if mean_counter==1
keep round regioncode gcounter dlns_mean dlnp unit_all unit *frac_hh* rain g2dum* g3dum* weight_counter
save `i'_savedmeans, replace
}




use grains_savedmeans, clear
label variable rain "rice"
label variable g2dum2_rain "chira"
label variable g2dum3_rain "khoi, lawa"
label variable g2dum4_rain "muri"
label variable g2dum5_rain "other rice products"
label variable g2dum6_rain "wheat"
label variable g2dum7_rain "maida"
label variable g2dum8_rain "suji, rawa"
label variable g2dum9_rain  "sewai noodles"
label variable g2dum10_rain "bread (bakery)"
label variable g2dum11_rain "other wheat products"
label variable g2dum12_rain "jowar"
label variable g2dum13_rain "bajra"
label variable g2dum14_rain "maize"
label variable g2dum15_rain "barley"
label variable g2dum16_rain "small millets"
label variable g2dum17_rain "ragi"
label variable g2dum18_rain "other cereals/cereal substitues"

save grains_savedmeans, replace

use veg_savedmeans, clear
label variable rain "potato"
label variable g2dum2_rain "onion"
label variable g2dum3_rain "radish"
label variable g2dum4_rain "carrot"
label variable g2dum5_rain "turnip"
label variable g2dum6_rain "beet"
label variable g2dum7_rain "sweet potato"
label variable g2dum8_rain "arum"
label variable g2dum9_rain "pumpkin"
label variable g2dum10_rain "gourd"
label variable g2dum11_rain "bitter gourd"
label variable g2dum12_rain "cucumber"
label variable g2dum13_rain "parwal, patal"
label variable g2dum14_rain "jhinga, torai"
label variable g2dum15_rain "snake gourd"
label variable g2dum16_rain "cauliflower"
label variable g2dum17_rain "cabbage"
label variable g2dum18_rain "brinjal"
label variable g2dum19_rain "lady's finger"
label variable g2dum20_rain "palak/othger leafy vegetables"
label variable g2dum21_rain "french beans, barbati"
label variable g2dum22_rain "tomato"
label variable g2dum23_rain "peas"
label variable g2dum24_rain "chillis: green"
label variable g2dum25_rain "capsicum"
label variable g2dum26_rain "plantain: green"
label variable g2dum27_rain "jackgruit: green"
label variable g2dum28_rain "lemon"
label variable g2dum29_rain "other vegetables"
save veg_savedmeans, replace



foreach i in grains veg{

use `i'_savedmeans, clear



**Hausman-Nevo instrument
bysort round gcounter: egen dlnp_other=mean(dlnp)
bysort round gcounter: gen counter=_N
gen state=floor(regioncode/10)
bysort round gcounter state: egen dlnp_state=mean(dlnp)
bysort round gcounter state: gen scounter=_N

replace dlnp_other=((dlnp_other*counter)-(dlnp_state*scounter))/(counter-scounter)


egen unit_rg=group(round gcounter)
egen unit_rr=group(round regioncode)

drop bfrac_hh

local estimator cue
local cutoff=0.66
local clus cluster unit_rr

local weight fweight
local controls dfrac_hh frac_hh
local high more_ru


local sample if (round==38 | round==66)



****Table 5: Elasticity estimates****

*OLS
areg dlns_mean dlnp `othendog' `othexog' `sample' [`weight'=weight_counter], absorb(unit_all) vce(`clus')
summ dlns_mean dlnp `othendog' `othexog' `sample'
cap: outreg2 using rsigma`i'_select`high'.tex, tex addstat(Adj R-squared,e(r2_a),"clusters",e(N_clust))  keep(dlnp  `othendog') replace ctitle(OLS)

areg dlns_mean dlnp `othendog' `othexog' dfrac_hh  `sample' [`weight'=weight_counter], absorb(unit_all) vce(`clus')
cap: outreg2 using rsigma`i'_select`high'.tex, tex addstat(Adj R-squared,e(r2_a),"clusters",e(N_clust))  keep(dlnp  dfrac_hh `othendog') append ctitle(OLS+dfrac)

areg dlns_mean dlnp `othendog' `othexog' `sample' & frac_hh>`cutoff' [`weight'=weight_counter], absorb(unit_all) vce(`clus')
cap: outreg2 using rsigma`i'_select`high'.tex, tex addstat(Adj R-squared,e(r2_a),"clusters",e(N_clust))  keep(dlnp  `othendog') append ctitle(OLS `cutoff')

*IV
**first-stage:
areg dlnp rain g2dum* `sample' [`weight'=weight_counter], absorb(unit_all) vce(`clus') 
cap: outreg2 using rsigma`i'_first.tex, tex label addstat(Adj R-squared,e(r2_a),"clusters",e(N_clust),"Joint F-stat",e(F))   keep(rain g2dum*) replace ctitle(rain)


reghdfe dlns_mean `othexog' (dlnp `othendog' = rain g2dum* `othiv') `sample' [`weight'=weight_counter], absorb(unit_all) vce(`clus') estimator(`estimator')
cap: outreg2 using rsigma`i'_select`high'.tex, tex addstat(Adj R-squared,e(r2_a),"clusters",e(N_clust))   keep(dlnp  `othendog') append ctitle(IV)


reghdfe dlns_mean `othexog' dfrac_hh (dlnp `othendog' = rain g2dum* `othiv') `sample' [`weight'=weight_counter], absorb(unit_all) vce(`clus') estimator(`estimator')
cap: outreg2 using rsigma`i'_select`high'.tex, tex addstat(Adj R-squared,e(r2_a),"clusters",e(N_clust))  keep(dlnp  dfrac_hh `othendog') append ctitle(IV+dfrac)


reghdfe dlns_mean `othexog' (dlnp `othendog'  = rain g2dum* `othiv') `sample' & frac_hh>`cutoff' [`weight'=weight_counter], absorb(unit_all) vce(`clus') estimator(`estimator')
cap: outreg2 using rsigma`i'_select`high'.tex, tex addstat(Adj R-squared,e(r2_a),"clusters",e(N_clust))  keep(dlnp `controls' `othendog') append ctitle(IV `cutoff')



*IV2
*First stage
areg dlnp dlnp_other  `sample' [`weight'=weight_counter], absorb(unit_all) vce(`clus') 
cap: outreg2 using rsigma`i'_first.tex, tex label addstat(Adj R-squared,e(r2_a),"clusters",e(N_clust),"Joint F-stat",e(F))   keep(dlnp_other) append ctitle(other prices)

reghdfe dlns_mean `othexog' (dlnp `othendog' = dlnp_other `othiv2') `sample' [`weight'=weight_counter], absorb(unit_all) vce(`clus') estimator(`estimator')
cap: outreg2 using rsigma`i'_select`high'.tex, tex addstat(Adj R-squared,e(r2_a),"clusters",e(N_clust))   keep(dlnp  `othendog') append ctitle(IV2)

reghdfe dlns_mean `othexog' dfrac_hh (dlnp  `othendog' =  dlnp_other `othiv2') `sample' [`weight'=weight_counter], absorb(unit_all) vce(`clus') estimator(`estimator')
cap: outreg2 using rsigma`i'_select`high'.tex, tex addstat(Adj R-squared,e(r2_a),"clusters",e(N_clust))  keep(dlnp dfrac_hh `othendog') append ctitle(IV2+dfrac)

reghdfe dlns_mean `othexog' (dlnp `othendog' =  dlnp_other `othiv2') `sample' & frac_hh> `cutoff' [`weight'=weight_counter], absorb(unit_all) vce(`clus') estimator(`estimator')
cap: outreg2 using rsigma`i'_select`high'.tex, tex addstat(Adj R-squared,e(r2_a),"clusters",e(N_clust))  keep(dlnp `controls' `othendog') append ctitle(IV2  `cutoff')

}







****Appendix: Testing for equality above/below median real per capita expenditure [Table A.8]

















*drop if mpce>=medpce

foreach i in grains veg{

local round

use temp`i'38, clear
append using temp`i'66


gen year=1983 if round==38
replace year=2009 if round==66

**merge in rainfall
merge m:1 regioncode year using "D:\Dropbox\indiadistrict\region_rain"
gen rain=act_rain_gca-norm_rain_gca

**region drop for rainfall analysis
keep if rain~=.
keep if _merge==3
drop if regioncode==.

*rural only
keep if sector==1
drop if reglat==. | reglon==.


bysort regioncode round: gen counta=_n
replace counta=0 if counta~=1
by regioncode: egen totcounta=total(counta)
drop counta


*threshold for share of housheolds purchasing
drop rnumhh
bysort round regioncode gcounter: egen frac_hh=total(hpresence)
bysort round regioncode gcounter: gen rnumhh=_N
replace frac_hh=frac_hh/rnumhh


drop rx rtotx
bysort round regioncode gcounter: egen rx=total(x)
bysort round regioncode: egen rtotx=total(x)
gen rshare=rx/rtotx

**drop households with missing data/implausible values
egen hh=group(round regioncode fsu hhid)

drop if x`i'==0
drop if hhsize>14
drop if mpce>10000 | mpce<10
drop if x==0
gen loghhn=log(n`i')

**only keep households with at least 2 varieties
drop if loghhn==0 | loghhn==.


compress



gen share=x/x`i'

replace mpce=mpce/7.256 if round==66
egen medpce=median(mpce)
drop if mpce<medpce

keep n`i'tot round* regioncode rn gcounter hpresence share rmeduv x`i' hhsize rain hh frac_hh reglat reglon mpce totcounta



egen unit_all=group(gcounter regioncode)
bysort round regioncode gcounter: gen rcounter=_n


levelsof n`i'tot, local(nmax)



gen state=floor(regioncode/10)
gen lns=log(share)


**make base the most popular variety in each region,

bysort regioncode gcounter: egen totx=max(frac_hh)
by regioncode: egen maxx=max(totx)
gen gmax=gcounter if totx==maxx
by regioncode: egen gbase=min(gmax)
summ gbase


gen blns=lns if gcounter==gbase
bysort hh: egen base_lns=min(blns)
gen brmeduv=rmeduv if gcounter==gbase
bysort hh: egen base_rmeduv=min(brmeduv)
gen bfrac_hh=frac_hh if gcounter==gbase
bysort hh: egen base_frac_hh=min(bfrac_hh)


gen dlns=lns-base_lns if lns~=0 & base_lns~=0 & gcounter~=gbase
gen dlnp=log(rmeduv)-log(base_rmeduv) if lns~=0 & base_lns~=0 & gcounter~=gbase
gen dfrac_hh=frac_hh-base_frac_hh if lns~=0 & base_lns~=0 & gcounter~=gbase




cap: drop base_* 
cap: drop blns 
cap: drop brmeduv 
cap: drop gmax
cap: drop maxx
cap: drop totx
cap: drop blnrs 

gen lnxpc=log(x`i'/(hhsize))
gen lnhhsize=log(hhsize)

foreach j in rmeduv{
replace `j'=log(`j')
}

gen prmeduv=rmeduv



foreach j in 66{
gen round`j'=0
replace round`j'=1 if round==`j'
}


egen region=group(regioncode)
egen regionmax=max(region)



sort round regioncode 
local regmax=regionmax[1]
forvalues j=1(1)`regmax'{
gen raintemp=rain if region==`j'
gen lat=reglat if region==`j'
gen lon=reglon if region==`j'
egen tlat=min(lat)
egen tlon=min(lon)
geodist reglat reglon tlat tlon, gen(dist`j') sphere
by round: egen rain_`j'=min(raintemp)
by round regioncode: egen distrain_`j'=min(rain_`j'/(1+log(1+dist`j')))
drop raintemp dist`j' lat lon tlat tlon rain_`j'
}



gen logrn=log(rn)


forvalues j=1(1)`nmax'{
gen pres=1 if gcounter==`j' & hpresence==1
bysort hh: egen pres`j'=min(pres)
drop pres
replace pres`j'=0 if pres`j'==.

gen price=prmeduv if gcounter==`j'
bysort round regioncode: egen price`j'=min(price)
bysort round state: egen price`j'_s=max(price)
bysort round: egen price`j'_ai=max(price)
replace price`j'=price`j'_s if price`j'==.
replace price`j'=price`j'_ai if price`j'==.
drop price price`j'_s price`j'_ai
}






egen unit_ind=group(regioncode)

xtset unit_all


cap: drop totcount
gen count=0
replace count=1 if lns~=.
bysort unit_all count round: gen r2counter=_n
replace r2counter=0 if r2counter~=1
bysort unit_all count: egen totcount=total(r2counter)


gen select=0
replace select=1 if dlns~=.


gen dfrac_hh2=dfrac_hh^2
gen dfrac_hh3=dfrac_hh^3
gen frac_hh2=frac_hh^2
gen frac_hh3=frac_hh^3


cap: drop zzz*
cap: drop gdum*
cap: drop g2dum*


drop if dlns==.


**construct a.good * rain dummies  AND   b. state * good *  rain dummies

egen zzz=group(gcounter)  

egen zzzmax=max(zzz)
local zzmax=zzzmax[1]
forvalues j=2(1)`zzmax'{
gen gdum`j'=0
replace gdum`j'=1 if zzz==`j'
gen g2dum`j'_rain=gdum`j'*rain
}
drop zzz zzzmax gdum*

egen zzz=group(gcounter state)  
egen zzzmax=max(zzz)
local zzmax=zzzmax[1]
forvalues j=2(1)`zzmax'{
gen gdum`j'=0
replace gdum`j'=1 if zzz==`j'
gen g3dum`j'_rain=gdum`j'*rain
}


bysort round regioncode gcounter: egen dlns_mean=mean(dlns)
by round regioncode gcounter: gen mean_counter=_n
by round regioncode gcounter: gen weight_counter=_N
egen unit=group(round regioncode gcounter)

keep if mean_counter==1
keep round regioncode gcounter dlns_mean dlnp unit_all unit *frac_hh* rain g2dum* g3dum* weight_counter
save `i'_savedmeans_high, replace
}





foreach i in grains veg{

local round

use temp`i'38, clear
append using temp`i'66


gen year=1983 if round==38
replace year=2009 if round==66

**merge in rainfall
merge m:1 regioncode year using "D:\Dropbox\indiadistrict\region_rain"
gen rain=act_rain_gca-norm_rain_gca

**region drop for rainfall analysis
keep if rain~=.
keep if _merge==3
drop if regioncode==.

*rural only
keep if sector==1
drop if reglat==. | reglon==.


bysort regioncode round: gen counta=_n
replace counta=0 if counta~=1
by regioncode: egen totcounta=total(counta)
drop counta


*threshold for share of housheolds purchasing
drop rnumhh
bysort round regioncode gcounter: egen frac_hh=total(hpresence)
bysort round regioncode gcounter: gen rnumhh=_N
replace frac_hh=frac_hh/rnumhh


drop rx rtotx
bysort round regioncode gcounter: egen rx=total(x)
bysort round regioncode: egen rtotx=total(x)
gen rshare=rx/rtotx

**drop households with missing data/implausible values
egen hh=group(round regioncode fsu hhid)

drop if x`i'==0
drop if hhsize>14
drop if mpce>10000 | mpce<10
drop if x==0
gen loghhn=log(n`i')

**only keep households with at least 2 varieties
drop if loghhn==0 | loghhn==.


compress



gen share=x/x`i'

replace mpce=mpce/7.256 if round==66
egen medpce=median(mpce)
drop if mpce>=medpce

keep n`i'tot round* regioncode rn gcounter hpresence share rmeduv x`i' hhsize rain hh frac_hh reglat reglon mpce totcounta



egen unit_all=group(gcounter regioncode)
bysort round regioncode gcounter: gen rcounter=_n


levelsof n`i'tot, local(nmax)



gen state=floor(regioncode/10)
gen lns=log(share)


**make base the most popular variety in each region,

bysort regioncode gcounter: egen totx=max(frac_hh)
by regioncode: egen maxx=max(totx)
gen gmax=gcounter if totx==maxx
by regioncode: egen gbase=min(gmax)
summ gbase


gen blns=lns if gcounter==gbase
bysort hh: egen base_lns=min(blns)
gen brmeduv=rmeduv if gcounter==gbase
bysort hh: egen base_rmeduv=min(brmeduv)
gen bfrac_hh=frac_hh if gcounter==gbase
bysort hh: egen base_frac_hh=min(bfrac_hh)


gen dlns=lns-base_lns if lns~=0 & base_lns~=0 & gcounter~=gbase
gen dlnp=log(rmeduv)-log(base_rmeduv) if lns~=0 & base_lns~=0 & gcounter~=gbase
gen dfrac_hh=frac_hh-base_frac_hh if lns~=0 & base_lns~=0 & gcounter~=gbase




cap: drop base_* 
cap: drop blns 
cap: drop brmeduv 
cap: drop gmax
cap: drop maxx
cap: drop totx
cap: drop blnrs 

gen lnxpc=log(x`i'/(hhsize))
gen lnhhsize=log(hhsize)

foreach j in rmeduv{
replace `j'=log(`j')
}

gen prmeduv=rmeduv



foreach j in 66{
gen round`j'=0
replace round`j'=1 if round==`j'
}


egen region=group(regioncode)
egen regionmax=max(region)



sort round regioncode 
local regmax=regionmax[1]
forvalues j=1(1)`regmax'{
gen raintemp=rain if region==`j'
gen lat=reglat if region==`j'
gen lon=reglon if region==`j'
egen tlat=min(lat)
egen tlon=min(lon)
geodist reglat reglon tlat tlon, gen(dist`j') sphere
by round: egen rain_`j'=min(raintemp)
by round regioncode: egen distrain_`j'=min(rain_`j'/(1+log(1+dist`j')))
drop raintemp dist`j' lat lon tlat tlon rain_`j'
}



gen logrn=log(rn)


forvalues j=1(1)`nmax'{
gen pres=1 if gcounter==`j' & hpresence==1
bysort hh: egen pres`j'=min(pres)
drop pres
replace pres`j'=0 if pres`j'==.

gen price=prmeduv if gcounter==`j'
bysort round regioncode: egen price`j'=min(price)
bysort round state: egen price`j'_s=max(price)
bysort round: egen price`j'_ai=max(price)
replace price`j'=price`j'_s if price`j'==.
replace price`j'=price`j'_ai if price`j'==.
drop price price`j'_s price`j'_ai
}






egen unit_ind=group(regioncode)

xtset unit_all


cap: drop totcount
gen count=0
replace count=1 if lns~=.
bysort unit_all count round: gen r2counter=_n
replace r2counter=0 if r2counter~=1
bysort unit_all count: egen totcount=total(r2counter)


gen select=0
replace select=1 if dlns~=.


gen dfrac_hh2=dfrac_hh^2
gen dfrac_hh3=dfrac_hh^3
gen frac_hh2=frac_hh^2
gen frac_hh3=frac_hh^3


cap: drop zzz*
cap: drop gdum*
cap: drop g2dum*


drop if dlns==.


**construct a.good * rain dummies  AND   b. state * good *  rain dummies

egen zzz=group(gcounter)  

egen zzzmax=max(zzz)
local zzmax=zzzmax[1]
forvalues j=2(1)`zzmax'{
gen gdum`j'=0
replace gdum`j'=1 if zzz==`j'
gen g2dum`j'_rain=gdum`j'*rain
}
drop zzz zzzmax gdum*

egen zzz=group(gcounter state)  
egen zzzmax=max(zzz)
local zzmax=zzzmax[1]
forvalues j=2(1)`zzmax'{
gen gdum`j'=0
replace gdum`j'=1 if zzz==`j'
gen g3dum`j'_rain=gdum`j'*rain
}


bysort round regioncode gcounter: egen dlns_mean=mean(dlns)
by round regioncode gcounter: gen mean_counter=_n
by round regioncode gcounter: gen weight_counter=_N
egen unit=group(round regioncode gcounter)

keep if mean_counter==1
keep round regioncode gcounter dlns_mean dlnp unit_all unit *frac_hh* rain g2dum* g3dum* weight_counter
save `i'_savedmeans_low, replace
}





local high _highlow

foreach i in grains{

use `i'_savedmeans_high, clear
gen high=1
append using `i'_savedmeans_low



**Hausman-Nevo instrument
bysort round gcounter: egen dlnp_other=mean(dlnp)
bysort round gcounter: gen counter=_N
gen state=floor(regioncode/10)
bysort round gcounter state: egen dlnp_state=mean(dlnp)
bysort round gcounter state: gen scounter=_N

replace dlnp_other=((dlnp_other*counter)-(dlnp_state*scounter))/(counter-scounter)


replace high=0 if high==.
gen high_dlnp=high*dlnp
gen high_rain=rain*high
foreach k of varlist g2dum2_rain-g3dum2_rain dlnp_other dfrac_hh{
gen high_`k'=high*`k'
}
drop high_g3dum2_rain

drop bfrac_hh


local cutoff=0.66
local clus cluster unit
local weight fweight
local controls dfrac_hh frac_hh 
local estimator cue


local othendog high_dlnp
local othexog high
local othiv high_rain high_g2dum*

local othiv2 high_dlnp_other


local sample if (round==38 | round==66) 



*OLS
areg dlns_mean dlnp `othendog' `othexog' `sample' [`weight'=weight_counter], absorb(unit_all) vce(`clus')
cap: outreg2 using rsigma`i'_select`high'.tex, tex addstat(Adj R-squared,e(r2_a),"clusters",e(N_clust))  keep(dlnp  `othendog') replace ctitle(OLS)

areg dlns_mean dlnp `othendog' `othexog' dfrac_hh high_dfrac_hh `sample' [`weight'=weight_counter], absorb(unit_all) vce(`clus')
cap: outreg2 using rsigma`i'_select`high'.tex, tex addstat(Adj R-squared,e(r2_a),"clusters",e(N_clust))  keep(dlnp  `othendog') append ctitle(OLS+dfrac)


*IV2
reghdfe dlns_mean `othexog' (dlnp `othendog' = dlnp_other `othiv2') `sample' [`weight'=weight_counter], absorb(unit_all) vce(`clus') estimator(`estimator')
cap: outreg2 using rsigma`i'_select`high'.tex, tex addstat(Adj R-squared,e(r2_a),"clusters",e(N_clust))   keep(dlnp  `othendog') append ctitle(IV2)

reghdfe dlns_mean `othexog' dfrac_hh high_dfrac_hh (dlnp `othendog' =  dlnp_other `othiv2') `sample' [`weight'=weight_counter], absorb(unit_all) vce(`clus') estimator(`estimator')
cap: outreg2 using rsigma`i'_select`high'.tex, tex addstat(Adj R-squared,e(r2_a),"clusters",e(N_clust))  keep(dlnp  `othendog') append ctitle(IV2+dfrac)
}




foreach i in veg{
use `i'_savedmeans_high, clear
gen high=1
append using `i'_savedmeans_low




**Hausman-Nevo instrument
bysort round gcounter: egen dlnp_other=mean(dlnp)
bysort round gcounter: gen counter=_N
gen state=floor(regioncode/10)
bysort round gcounter state: egen dlnp_state=mean(dlnp)
bysort round gcounter state: gen scounter=_N

replace dlnp_other=((dlnp_other*counter)-(dlnp_state*scounter))/(counter-scounter)


replace high=0 if high==.
gen high_dlnp=high*dlnp
gen high_rain=rain*high
foreach k of varlist g2dum2_rain-g3dum2_rain dlnp_other dfrac_hh{
gen high_`k'=high*`k'
}
drop high_g3dum2_rain

drop bfrac_hh


local cutoff=0.66
local clus cluster unit
local weight fweight
local controls dfrac_hh frac_hh 
local estimator cue

local othendog high_dlnp
local othexog high
local othiv high_rain high_g2dum*

local othiv2 high_dlnp_other


local sample if (round==38 | round==66) 

*OLS
areg dlns_mean dlnp `othendog' `othexog' `sample' [`weight'=weight_counter], absorb(unit_all) vce(`clus')
cap: outreg2 using rsigma`i'_select`high'.tex, tex addstat(Adj R-squared,e(r2_a),"clusters",e(N_clust))  keep(dlnp  `othendog') replace ctitle(OLS)

areg dlns_mean dlnp `othendog' `othexog' dfrac_hh high_dfrac_hh `sample' [`weight'=weight_counter], absorb(unit_all) vce(`clus')
cap: outreg2 using rsigma`i'_select`high'.tex, tex addstat(Adj R-squared,e(r2_a),"clusters",e(N_clust))  keep(dlnp  `othendog') append ctitle(OLS+dfrac)


*IV

reghdfe dlns_mean `othexog' (dlnp `othendog' = rain g2dum* `othiv') `sample' [`weight'=weight_counter], absorb(unit_all) vce(`clus') estimator(`estimator')
cap: outreg2 using rsigma`i'_select`high'.tex, tex addstat(Adj R-squared,e(r2_a),"clusters",e(N_clust))   keep(dlnp  `othendog') append ctitle(IV)

reghdfe dlns_mean dfrac_hh high_dfrac_hh `othexog' (dlnp `othendog' = rain g2dum* `othiv') `sample' [`weight'=weight_counter], absorb(unit_all) vce(`clus') estimator(`estimator')
cap: outreg2 using rsigma`i'_select`high'.tex, tex addstat(Adj R-squared,e(r2_a),"clusters",e(N_clust))  keep(dlnp  `othendog') append ctitle(IV+dfrac)
}





